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ABSTRACT 

We investigate through hydrodynamic simulations the destruction of ncwly-formcd dust grains by 
sputtering in the reverse shocks of supernova remnants. Using an ideahzed setup of a planar shock 
impacting a dense, spherical clump, we implant a population of Lagrangian particles into the clump to 
represent a distribution of dust grains in size and composition. We then post-process the simulation 
output to calculate the grain sputtering for a variety of species and size distributions. We explore the 
parameter space appropriate for this problem by altering the over-density of the ejecta clumps and 
the speed of the reverse shocks. Since radiative cooling could lower the temperature of the medium in 
which the dust is embedded and potentially protect the dust by slowing or halting grain sputtering, 
we study the effects of different cooling methods over the time scale of the simulations. In general, 
our results indicate that grains with radii less than 0.1 fj,m are sputtered to much smaller radii and 
often destroyed completely, while larger grains survive their interaction with the reverse shock. We 
also find that, for high ejecta densities, the percentage of dust that survives is strongly dependent on 
the relative velocity between the clump and the reverse shock, causing up to 50% more destruction for 
the highest velocity shocks. The fraction of dust destroyed varies widely across grain species, ranging 
from total destruction of AI2O3 grains to minimal destruction of Fe grains (only 20% destruction in 
the most extreme cases). C and Si02 grains show moderate to strong sputtering as well, with 38% 
and 80% mass loss. The survival rate of grains formed by early supernovae is crucial in determining 
whether or not they can act as the "dust factories" needed to explain high-redshift dust. 
Subject headings: hydrodynamics — supernova remnants — shock waves — dust, extinction 



1. INTRODUCTION 

The presence of dust grains in the interstellar medium 
(ISM) of galaxies can have a wide range of impacts on 
the evolution of the universe, both at present and at high 
redshift. As refractory elements become depleted from 
the gas phase and locked into grains, they form a con- 
siderable reservoir of coolants in the solid-state phase. 
The surfaces of these grains are also believed to act as 
formation sites for molecular hydrogen (II2), which not 
only serves as an additional coolant but can provide the 
initial conditions for a rich cloud chemistry. Grains are 
capable of reprocessing the radiation field of the ISM 
by providing a critical transfer mechanism for turning 
ultraviolet /optical starlight into far- infrared and submil- 
limeter emission seen in interstellar molecular clouds and 
high-redshift galaxies. The earliest dust grains may also 
have influenced the low-metallicity gas by modifying the 
thermodynamic conditions in the high-redshift universe 
and altering the star formation process and the initial 
mass function (IMF) for the earliest second -generation 
stars ( jOmukai et al.||2005[ park et al.||2008[ ). 

One paradox arises from observations of dust at high 
redshift. Several observations of da mped Lyg systems 



Gyr after the Big Bang, which is difficult to explain us- 
ing traditio nal models of grain formation in evolved low- 



mass stars (Dunne et al. 2003 Dwek et al. 20071 and 



transportation of dust into the iSM through stellar winds 
( |Whittet|[T992| . 

In an effort to find other sources for high-redshift dust, 
recent models have looked at the ejecta of core-collapse 
supernovae (SNe) as possible formation sites, which can 
occur on ti mescales shorter tha n the evolution of low- 
mass stars ( Dwek fc Scalo||T980 1 . Since the inner ejecta 
of core-collapse SINe are generally cold, dense, and metal- 
enriched, they are conducive to producing significant 
amounts of new dust grains. By utilizing theories of nu- 
cleatio n and grain g rowth made b y Kozasa & Hasegawa 
(11987'), two groups ( [Kozasa et al.||1989( |1991( |Todini ^ 
Ferrara 2001) have attempted to estimate the dust mass 
formed in expanding SNe ejecta of varying metallicities, 
with results ranging fr om ^0.1-0.3 Af^. Th ese i nitial re 



suits were extended by Nozawa et al. ( 2003 ) and Schnei 



der et al. (2004) to the case of zero-nietallicity SiNe with 



indicate the presenc e of dust grains ( Pettini et al.||1994 
Ledoux et al. 2002), and high-redshift [z > 6) quasars 
seen in the Sloan Digital Sky Survey (SDSS) along with 
daxies in the Hubble Deep Field show 1.2 mm ther- 



significant production of O, Si, and other heavy elements, 
as would be created by Population HI stars. 

Although SNe are appealing as potential "dust facto- 
ries" , current observations of nearby SNe appear to fall 
short of the required dust masse s (~0.2 M^^ or po ssi- 



mal emission from dust dSmail et al.||1997| Hughes et al. 


1998 


Bertoldi et al. 


20U3). 'I'hese observations imply 



a total dust mass of ~T0^ Mcd formed within the first 



bly as high as 1 M^^ per event; Dwek et al. 2007) by 
factors of 10-100. Sugerman et al. ~i 2006 ) produced a 
review of Type II SN 2003gd in whicn they used radia- 
tive transfer models to conclude that up to 0.02 A/0 of 
dust is present in t he ejecta, but this is contradicted by 
Meikle et al.|(|2007|) who found only 4 x 10~^Mo of newly 
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formed dust using mid-IR observations. Similar values 
were found in a young supernova remnant (SNR) l ocated 



in the Small Magella nic Cloud (SMC) studied by |Stan- 
who claimed to observe 10"'^ Mr 



( 2005 ) 



imirovic et al 

ot hot dust. Kotak et al. ( 2009 ) assert a direct detection 
of ejecta dust in I'ype 11-plateau SN 2004et, but again 
only estimate a few times 10~^ Mq. One promis ing ob- 
servat i on for SNe as dust factories carried out by |Dunne| 
et al. (2003) argued for the existence of ~3 Mq of dust 



in the direction of Cassiopeia A based on a sub millime 
ter detection , but it was quickly contradicted by |Kraus"e 



et al. ( 2004 ) who attributed this emission to interstellar 
dust contamed in an adjacent molecular cloud. In the 



(20061 



Large Magellanic Cloud (LMC), [Williams et al 
observed four SNRs using the 24 and 70 //m bands on 
Spitzer Space Telescope a nd inferred dust masses o f 0.01- 
0.1 Mq. More recently, |Rho et al.| ( |2008[ |2009| used 
Spitzer to study Cassiopeia A and young supernova rem- 
nant 1E0102-7219 to find 0.020-0.054 Mq and 0.014 M© 
of dust respectively. To continue the hunt for dust in 
Cassiopeia A, data from AKARI and the Balloon-borne 
Large A perture Submillimeter Telescope (BLAST) were 
used by jSibthorpe et al.| ( |2009[ ) to find 0.06 Mq of new 
"cold" ( i' ~ 35 K) dust, which when combined with pre- 
vious "warm" (T ~ lOOK) dust estimates, could make 
the dust-factory argument more plausible. 

Observational tests that seek to confirm the predic- 
tions made by models of SNe are difficult. It is neces- 
sary to disentangle far-IR backgrounds from the thermal 
dust emission and separate the newly formed ejecta dust 
from grains contained in circumstellar and interstellar 
material. One must also determine what fraction of the 
dust can be destroyed by the interaction between SNR 
reverse shocks and the ejecta. As a result of these dif- 
ficulties, some dust may be located in undetected, un- 
shocked ejecta. SN 1987A presents a particularly good 
case of the challenge of producing accurate dust mass es- 
timates, since it has been inferred that dust formation 
is almost certainly o ccurring, but only 10~^ Mq was de 



tected (Dwek 2006 Dwek & Arendt 2007). Observed 



changes in the optical and boiometric Muxes and asym- 
metric red and blue emission-lin e wings indica t e tha t 



more dust is most likely present (McCray 1993 2007) 
Furthermore, since SN 1987A is a young remnant, the 
reverse shock has not yet impacted the inner ejecta ma- 
terial — a location that is very likely to harbor additional 
dust. 

The question arises as to whether or not the fast- 
moving dust (Voj > 1000 km s~^) can survive its impact 
with the reverse shock, at which point it will be subject to 
thermal sputtering and grain-grain collisions. Since the 
relative speed between the shock and the ejecta clumps 
can be large, the temperature of the ejecta can be raised 
to a level that might completel y destroy these newly 



formed grains. Recent work by Nozawa et al. (20071 



investigates the interaction between reverse shocks and 
SN ejecta through one-dimensional hydrodynamic simu- 
lations. They conclude that 20-100% of the dust mass is 
destroyed, depending on the ambie nt gas density and the 
initial energy of the SN explosion. [Bianchi fc Schneider 



(2007) also used a combination of numerical and semi 
analytic methods to conclude that only a small fraction 
of the dust surives — as low as 7% in some of their mod- 



els. Further efforts were made by Nath et ah] ( 2008 ) who 
performed one-dimensional analytic calculations of the 
reverse shock impacting ejecta material that followed a 
power-law density profile. They found that only 1-20% 
of the dust mass in the form of silicates and graphites 
was sputtered, depending on the steepness of the profile 
and the grain size distribution. 

While the above studies have treated the ejecta ma- 
terial primarily as a continuous medium, observations 
of SNR ejecta indica te that in reality the material is 
knotted a nd clumpy (iFesen et al.||2006 Hammell & Fe 



sen 



2008[ ). This lends itself well to The classic "cloud- 
crushing problem ( jWoodward^^l976| ) in which an over- 
dense clump/cloud ot material is impacted by a shock 
wave. Such a problem has been studied using numer- 
ical hydrodynamic simulations under in the context of 
an ISM cloud getti ng hit by the forward shock of a su- 
pernova explos ion (Bedogn i fc Woodward|1990| Stone fc 
Norman] [T992ljK lcin et aT VM \ [Mac Low et al.||li)Ml 
Orlando et al.||2005i |Patnaude fc Fesen||2005 Nakamura 
et al.||2006| ). 



In this paper, we investigate the effects of reverse 
shocks on SN ejecta by performing three-dimensional 
cloud-crushing simulations and focusing on the case of 
an individual ejecta clump encountering a reverse shock. 
In Section [2j we describe the code that we have used 
to carry out these simulations and our particular simu- 
lation set up. In Section |3j we discuss the simulations 
examined in this paper, which include comparisons to 
previous work and runs that serve as our means of com- 
puting reverse-shock dust destruction. In Section |4j we 
highlight our methods for calculating dust destruction. 
In Section [5] we present the results of calculations, and 
wc conclucle in Section |6] with a summary and discussion 
of future work. 

2. METHODOLOGY 
2.1. Enzo 

For this study, we use the cosmological, Eulerian adap- 
tive mesh re finement (AMR), hydrodynamics -|- N-body 
code, Enzo ([Bryan fc Norman 1997 Norman fc Bryan 
19991 10 'Shea et al.|2005|). Owing to the idealized nature 



of our problem, we do not use any of the cosmological 
or gravity-solving components of the code. To handle 
the fluid dynamics in our work, we use En zo's picccwise 
parabolic method (PPM; [Colella fc Woodward,1984| ) for 
solving the hydrodynamic equations. l^'Gr more details 
on the inner-workings of Enzo, refer to the cited articles. 

The AMR capabilities of Enzo are particularly useful 
to resolve a range of physical scales or to ensure that only 
the regions of particular concern remain well resolved, so 
that we do not waste computation time on unimportant 
or uninteresting portions of the computational domain. 
The refinement algorithm works by taking any cells that 
meet a user-specified refinement criterion and spawning a 
new grid of cells with twice the resolution. For example, 
a three-dimensional cell that satisfies any of the desired 
refinement criteria will be split into eight cells with cell 
edges that are now half as long as the parent cell from 
which they were spawned. Once this cell is refined to a 
higher level, all of the new child cells are then checked 
to see if any still meet the refinement parameters. If so, 
they will be refined to a higher level in the same way, be- 
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coming parents to new children. Higher and higher levels 
of refinement can be created to theoretically infinite res- 
olution. However, since large amounts of refinement can 
quickly become computationally expensive, a maximum 
level is often specified when running to the code to halt 
refinement at the highest resolution that can be afforded 
by the simulator. Enzo comes preloaded with the ability 
to refine based on a variety of parameters: matter over- 
density; gradients in pressure, energy, or density; Jeans 
length; cooling time; and shocks. It is also possible to de- 
fine new refinement criteria that are simulation specific. 
We have created such a criterion for our cloud-crushing 
simulations that refines on gas that is initially contained 
within the cloud. The impleme ntat ion of that refinement 
method is discussed in Section [321 

In addition to the standard hydrodynamic capabilities 
of Enzo, wc utilize "tracer particles" in all of our simula- 
tions. The particles are Lagrangian in nature and simply 
advect with the flow present in the grids. In the current 
format, the tracer particles do not alter the state of their 
surrounding medium and serve only as a means of track- 
ing the conditions (i.e., density and temperature) of the 
fluid they are embedded in. However, the machinery for 
allowing the tracer particles to produce feedback, such as 
metal injection or dust cooling, exists and will be utilized 
in our future papers. Currently, the tracer particles are 
particularly useful in representing the presence of dust in 
our simul ation s, and their use as such will be addressed 
in Section |4]TJ 

For a subset of our simulations, we also include the 
effects of radiative cooli ng. Currently, we use t he ana- 
lytic formula defined by Sarazin & White (1987), which 
approximates the cooling rate as a function of tempera- 
ture for Z = 0.5Zq gas in ionization equilibrium. While 
we expect the ejecta material to be metal enriched, with 
relative abundances that likely deviate from solar values, 
for the purposes of this paper we sought only to inves- 
tigate the effect of the presence of metals (cooling on) 
versus a complete lack of metals (cooling off). Future 
work will include more variable abundances to allow for 
the enrichment of the ejecta by metals such as C, O, Si, 
Mg, Mn, and Fe which have been shown to have signifi- 
cant impacts on the cooling rates ( Santoro fc Shull|2006 1 
and are key ingredients in most dust grains. When imple- 
menting these adjustable abun dances and altered c ooling 
rates, we plan to use a Cloudy ( Ferland e t al.|1998 ) cool- 
ing method similar to Smith et al. (2008( ). 



2.2. Simulation Setup 

To address the proposed problem of a supernova rem- 
nant reverse shock impacting a clump of newly formed 
ejecta, we added to the suite of pre-packaged Enzo simu- 
lations. In this particular problem type, we seek to sim- 
ulate the previously studied "cloud-crushing" problem in 
which a shock is driven into a cloud of over-dense ma- 
terial. Although previous work focused on the forward 
shock of a supernova propagating into the ISM, the pro- 
cess is analogous to that of a reverse shock penetrating 
the clumpy supernova ejecta. 

The basic nature of our new problem type is to initial- 
ize a spherical cloud with constant overdensity, x, em- 
bedded in an ambient medium of density pm such that 
Pcioud = XPm- We then send a shock wave with Mach 
number, A4 = v/cs where Cs is the adiabatic sound speed 



with 7 = 5/3, in the cartesian z-direction toward the 
cloud. Additional user-defined parameters required to 
initialize the simulations include the radius of the cloud, 
''cloud, and the temperature of the cloud, Tdoud- AH of 
the remaining quantities needed to start the simulation 
are derived from the above mentioned user-specified pa- 
rameters. These derived values include temperature of 
the ambient medium, Tm, velocity of the shock, Wshock, 
post-shock density, Pshock, and post-shock temperature, 
Tshock- We set Tm so that the cloud remains in pres- 
sure equilibrium, while the shock-related values are cal- 
culated using the Rankine-Hugoniot jump conditions, 
which when defined in the frame of the shock front, are 
given by: 



(1) 



P2 




' (7 + l)M2 ■ 




Pi 




_{j-l)M^ + 2 




Vl 






(7+1)^(2 ■ 




V2 




(7- 1)M2 +2 












1) 




.P2 




[ (7+1) 



(2) 



(3) 



Here, the subscripts 1 and 2 refer to the upstream and 
downstream quantities of the flow with respect to the 
shock front and Ai — vi/cg- See Figure jl^ for a 
schematic of the flow. From these jump conditions, we 
shift from the rest frame of the shock front to the rest 
frame of the ambient medium to define the quantities 
used in our simulation. Under this transformation, the 
upstream and downstream values for density and tem- 
perature become the medium and shock quantities, re- 
spectively. However, the velocity of the inflowing ma- 
terial is defined as the difference between upstream and 
downstream values, Wshock = vi — V2- This results in the 
following set of equations to define the initial conditions 
for the simulation: 



Pshock — Pn 



t shock 



(7 + 1)7^2 ' 



277W2- (7- 1) 



Pshock 



(7+1) 



■1 



Pshock / 



(4) 



(5) 



(6) 



where the meanings of A4 and Cg remain unchanged and 
the fiow is now represented by Figure [T|d. 

While this formulation of the cloud-crushing problem 
is sufficient to study the general properties of the shock- 
cloud interaction and compare to previous works, we 
made provisions for additional complexities worthy of in- 
vestigation. In one avenue of investigation, we allowed 
the cloud to take on a variety of density profiles. While 
the simplest implementation is to define the cloud as a 
sphere of uniform density, which we explore in some of 
our simulations, we expect that the boundary between 
the ejecta clump and the ambient medium shoul d be 
smoothly varying, ana logo us to the ISM clouds of Pat- 
naude & Fesen (2005) and Nakamura et al. (2006). In 
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Fig. 1. — (a) Jump condition quantities in the frame of the shock front, represented by the vertical black line, (b) Depiction of cloud- 
crushing setup in the frame of the ambient medium. The medium and the cloud are both at rest. 



addition, a "soft" cloud boundary reduces undesirable 
numerical effects commonly produced by a hard edge. 
Therefore, we define the cloud to consist of a smoothly 
varying envelope surrounding a uniform core for a frac- 
tion of our simulations. Specifically, we allow for two 
types of cloud envelopes — one given by a 1/r power- 
law fall- off. 



Pc 



Pif) = { Pcorc + [Pm - Pc 



:0<r <rr 



-l/^cor. 



(7) 



7 ^corc <^ ^ ^cloud 

and another by a gaussian fall-off, 

Pcorc 1 ^ ^ ^ ^core 



p{r) = < 



Pc 



cloud~^'corc , 



(8) 



< r < r, 



loud 



where rcore is the inner radius at which the cloud is as- 
sumed to have a uniform density and Pcoro = XPm is the 
density of that core. In the case of a uniform cloud Pcorc 
and Pcioud are interchangeable as indicated in Figure [1]d. 

For all of the simulations presented in this paper, we 
have treated the shock as a constant infiow of material. 
This is representative of the scenario in which the width 
of the shock is much large than the radius of the cloud 
(effectively infinite) . For the situation in which the shock 
is of finite width and comparable in size to that of the 
cloud, it would be necessary to send a time-limited pulse 
of post-shock material. Such methods are being investi- 
gated for future work. 

2.3. Cloud Tracking 

Finally, we implement two additional features that al- 
low us to track the evolution of the ejecta clump. First, 
we populate all of our clouds with a uniform density of 
tracer particles to represent dust contained in the ejecta 
material. This allows us to track the movement of the 



dust while the cloud is crushed, as well as the nature of 
the medium in which the dust is embedded. We can then 
post-process the history of each tracer particle to calcu- 
late dust destruction by assigning each particle some dis- 
tribution of dust-grain properties. Second, we add a sim- 
ple scalar field that is advected with the flow in the same 
manner as density, but does not alter the state of the sys- 
tem. We use this field to represent "cloud material", ^, 
by initializing the values to be non-negligible only in cells 
that fall within the radius of the cloud. This field is mo- 
tivated purely by computational necessity, and all cloud 
material values are dimensionless and devoid of physical 
meaning. To avoid potential round-off errors in the nu- 
merics, we set cells that are outside the cloud radius to 
be small but non-zero (10^^''). The non-negligible val- 
ues are set to be dependent on the density of the cloud 
such that a cell with initial density pcorc is assigned a 
cloud material value of ^ = 1. This means that a uni- 
form density cloud would be filled with I's while a cloud 
with a power-law or gaussian envelope would have I's in 
the core and then a fall-off in cloud material that scales 
directly with the density profile. This field proves to be 
particularly useful as a refinement criterion for the AMR 
component of Enzo to ensure that the hydrodynamic evo- 
lution of th e cloud is tracked with sufficient resolution. 
Section 13.21 discusses its uses as a refinement criterion. 

3. THE SIMULATIONS 

3.1. Dimensionality and Comparison to Previous Works 

Before addressing the ejecta-reverse shock interaction, 
we present simulations in which we demonstrate Enzo's 
ability to successfully handle the cloud-crushing prob- 
lem. We also compare it to previous studies that inves- 
tigated the interactions between the forward shock and 
the ISM. In this simulation, we construct a computa- 
tional volume that consists of 512^ cells and a physical 
size of 2.47 x 10^^ cm (0.8 pc) on each side. We de- 
fine the cloud of over-dense material to have a radius 
''cloud = 3.09 X 10^'' cm (0.1 pc) and an over-density fac- 
tor X = 10 in a higher resolution nested grid such that 
there are 128 cells per cloud radius, a physical resolution 
of 2.41 X 10^^ cm (^161 AU) pe r ce ll edge. This choice is 



motivated by Klein et al. ( 1994 ) and Orlando et al. ( 2005 1 
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who found that this provides sufficient resolution to cap- 
ture the physical properties of the problem. We then use 
inflow boundary conditions to send material with the ap- 
propriate shock values to produce a = 10 shock wave 
in the direction of the cloud. We allow the simulation to 
run for a time t = 3icc, where tec = Y^^^?-nioi.ri /«shnnk i s 
the cloud-crushing time as defined by [Klein et al. (1994), 
at which point the cloud flows out of the domain! We 
find that this set of input parameters compares closely t o 
simulation s carried out by Bedo gni fc Woodwa rd (1990) , 
Stone fc Norman] ( 119921) , Klein ct al.| ( |1994[ )7ahd Orlando 
et al. ( j2005[ ) ana thereiorc creates a simulation that can 
be compared to those previous works. We maintain the 
2.41 X 10^^ cm per cell edge resolution over the computa- 
tional domain that encases the cloud for the full duration 
of the simulation. 

In order to explore the effects of dimensionality on this 
problem, we investigate the evolution of this particular 
set up by starting first with a one-dimensional (ID) sim- 
ulation and then progressing to higher dimensions. The 
ID case corresponds to a scenario in which the inflowing 
shock impacts an effective "wall" of cloud material. This 
situation lends itself well to tracking the various trans- 
mitted and reflected shocks present in the problem as 
well as for identifying these features in higher dimension 
simulations. The first example of such shock propagation 
occurs when the infiowing shock impacts the cloud wall 
and transmits a slower shock through the cloud material 
and a bow shock is reflected back towards the inflow. 
When the first transmitted shock crosses the back end 
of the cloud wall, another shock is transmitted into the 
ambient medium, while a rarefaction wave travels in the 
direction opposite the inflow through the newly-shocked 
cloud material. At this point, additional ringing occurs 
as waves bounce back and forth inside the cloud wall. 
See Figure [2] for density and temperature profiles of the 
ID simulation at various points in time which display the 
above mentioned features. 

For the two-dimensional (2D) case, we find many of the 
same features as in the ID case, with the exception that 
the infiowing shock is capable of wrapping around and re- 
converging behind the cloud before the shock transmitted 
into the cloud traverses the full extent of the cloud. This 
is because the cloud no longer has the wall-like presence 
it had in the ID case, and is now a circular object sur- 
rounded by ambient medium. As a result, an additional 
reverse shock is driven through the cloud opposite the 
flow as the inflowing shock envelops the cloud and exerts 
added pressure on all sides. The second dimension also 
allows for vortex creation, as material is stripped off the 
outer edges of the cloud and pushed into the evacuated 
region directly behind the cloud by the high pressure in- 
flow material that surrounds it. In our 2D realization of 
this pro blem, we find qualitatively sim ilar results to the 
work of Bedogni & Woodward ( 1990 ) — further confir- 
mation that Fnzo provides adequate machinery to study 
this problem. Figure [3] shows the results of this 2D run 
in density, temperature, and enstrophy, the square of the 
magnitude of the vorticity vector, |wp = |V x vp. 

In Figure |4j we see that the time evolution of the pro- 
jected density for the three-dimensional (3D) simulation 
also compares favorably with previous works, particu- 
larly when we track the development of various features 
as a function of time. As in the ID and 2D cases, when 



the inflowing shock first impacts the cloud, we see a 
slower shock transmitted into the cloud and a reflected 
bow shock that propagates into the post-shock medium 
while the cloud becomes laterally fiattened. By time 
t — 0.75icci the shock has flowed completely past the 
cloud and converged directly behind the cloud on the 
X-axis. At this point, the cloud is bathed in the high- 
pressure, post-shock medium and is compressed in all 
directions; an additional shock is driven back into the 
cloud opposite the direction of the inflow. Eventually, 
around t = 1.2tcc, the cloud begins to re-expand into 
the post shock medium owing to the compression caused 
when the shock first passed over it. Finally, for t > 2.0tcc, 
the cloud is dominated by hydrodynamic instabilities and 
takes on a highly non-uniform, fragmented, and filamen- 
tary structure as it becomes comoving with the inflowing 
material. The evolution of vorticity in the 3D simula- 
tion becomes even more complex than the 2D case, as 
the third dimension allows for the stretching of vortex 
tubes and leads ultimately to a greater level of symmetry- 
breaking turbuleiice. A l l of these features wer e seen by 
Stone fc Norman] ( |1992[ ), [Klein et al.| ( |1994[ ), and ]Orlando 
et al. "( |?005 ) at comparable moments m time. Figure^ 



also shows the same evolution in projected temperature 
To produce a more direct comparison of the three cases 
of varying dimensionality, we found it useful to look at ID 
proflles of density and temperature as a function of time. 
For the ID data, these profiles are plots of the data as in 
Figure [2] while for the 2D and 3D simulations, they are 
found by tracing a path through the center of the compu- 
tational domain that runs parallel to the a;-axis. Figure[5] 
presents these profiles for select time steps. Upon inspec- 
tion, the 2D and 3D cases are considerably more complex 
than the ID case, but they still show similar evolution, 
especially in the time period before the inflowing shock 
wraps around the cloud in the 2D and 3D cases. We 
also point out that the cloud material moves across the 
computational domain more rapidly for lower dimension 
simulations. This can be explained by the fact that the 
momentum imparted by the planar shock has fewer de- 
grees of freedom to transport into. For the ID case, all 
of the material is forced to move in the x-direction by 
the inflow, but in the 2D and 3D cases it can be diffused 
into the y- and z-directions. 

3.2. Reverse Shocks and Ejecta Clumps 

While the simulations we use to study the interaction 
between a clump of ejecta and a supernova remnant re- 
verse shock are qualitatively similar to the case of a for- 
ward shock penetrating the ISM, the values of the input 
parameters are considerably different. In the case of the 
reverse shock impacting newly formed ejecta, the relative 
velocity between the cloud and the shock wave is gener- 
ally an order of magnitude higher (^lO'^ km s~^) than 
the shock-ISM case (~10^ km s~^), and the overdensity 
of the ejecta clump is on the order of x ^ 1000 versus 
the X ~ 10 values assumed for ISM clouds. 

In addition to the modifications of the input parame- 
ters, we also made use of the variable cloud structure to 
explore what role, if any, the density profile of the cloud 
played in determining the amount of dust that would 
be sputtered as the ejecta are shocked. Using the full 
array of configurable options in our cloud-crushing prob- 
lem type (shock width, cloud structure, radiative cooling, 
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Fig. 2. — Density (top) and temperature (bottom) profiles for the ID comparison simulation. The left column corresponds to t = 0.52tcc, 
and the right column corresponds to t = 1.40tcc, where 4 = is when the shock first impinges upon the cloud edge. The labels BS, TSl, 
TS2, and RW mark the locations of the bow shock, first transmitted shock, second transmitted shock, and rarefaction wave respectively. 



cloud overdensity, shock velocity), we ran a grid of sim- 
ulations to determine the parameters that strongly in- 
fluence the dust destruction induced by the shock-cloud 
interaction. Specifically, we investigate various combina- 
tions of power-law and gaussian envelopes; cloud over- 
densities of X = 100 and 1000; shock velocities of 1000, 
3000, and 500 km s~^; and no coo ling and the ana- 



lytic cooling of Sarazin & White ( 1987| . For a list of the 
simulations and their corresponding parameters refer to 
Table HI 

All of the simulations have root-grid dimensions of 
512x256x256, which provides a physical resolution of 
1.25 X 10^^ cm (~84 AU) per cell edge. At this level, the 
radius of the cloud spans 8 cells for a physical radius of 
10^^ cm. In addition to the base-level resolution, we in- 
clude a problem-specific refinement criterion that allows 
for two additional levels of resolution. Cell refinement is 
controlled using the previously mentioned cloud material 
field and flags cells for refinement based on a dimension- 
less cloud material "mass", = ^Vcdi, where Vceii is 
the volume of the cell. To determine whether or not a 
cell should be flagged it must have a cloud material mass 
greater than mfiag = Timin'''") where mj^in is the mini- 
mum mass for refinement, r is the refinement factor, I is 
the current refinement level of the cell, and a is a tun- 
able exponent used to control whether the refinement 
acts in a Lagrangian (a — 0), super-Lagrangian (a < 0), 



or sub-Lagrangian (a > 0) manner. For all simulations 
presented in this paper, we use TOmin = 10^^, r = 2, 
and a = 0, where m^in is a dimensionless number. By 
allowing the grid to refine to two levels beyond the root 
grid, we achieve a physical resolution of 3.125 x lO^'* cm 
(~21 AU) on the highest level. Although these reverse- 
shock simulations do not have the same number of cells 
per cloud radius as presented in the comparison case, 
we find that this particular combination of initial root- 
grid resolution, cloud material refinement criterion, and 
two-level refinement is sufficient to track and study the 
properties of the interactions between the reverse shock 
and the ejecta cloud. Additionally, the final dust masses 
do not change significantly with increased resolution. 

We set the duration of our reverse-shock simula- 
tions based on the longest computationally manageable 
timescale for the parameters that produce the longest 
cloud-crushing time. For the above-mentioned parame- 
ter space, this corresponds to simulations with Tdoud = 
10^6 cm, X = 1000, and Wshock = 1000 km s'^, giv- 
ing tec ~ 100 yrs. With this constraint, we are able to 
run these simulations to a physical time of t 10'^ yrs 
~ lO^cc- Therefore, we run all simulations to a mini- 
mum of 10 cloud-crushing times which results in varied 
physical times. 

Owing to the nature of the inflowing shock and rele- 
vant time scales, we encountered a possible simulation 
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Fig. 3. — Density (top), temperature (middle), and enstrophy (bottom) for the 2D comparison simulation. The left column shows the 
state of the system at i = 0. Slice, the middle at t = 1.69tcc, and the right at t = 2.56tcc 



limitation in which the cloud material was often blown 
out of the back end of the computational domain by the 
shock medium. We solved this limitation and kept the 
majority of the cloud material in the simulation box for 
the desired duration by allowing the simulations to run 
up to a point in which the cloud is nearly comoving with 
the inflowing shock material and then shift into the in- 
flow frame. We execute this by pausing the calculation 
at a suitable moment in time, subtracting off the inflow 
velocity, Wghock: from the x- velocity of every cell, and re- 
suming the simulation. After the simulation is restarted, 
the cloud evolves mostly statically in terms of bulk mo- 
tion. Since the temperatures and densities remain unal- 
tered, the cloud's relative motions are identical to what 
they would have been in the initial rest frame. 

In studying the evolution of these "ejecta-crushing" 
simulations, we find that at early times the simula- 
tions agree with the previously discussed 3D compar- 



ison simulation, with the exception that the smoothly 
varying envelope subdues the interactions between the 
cloud boundary and the inflowing material. The lack of a 
sharp density discontinuity decreases the development of 
features perpendicular to the direction of the inflow and 
the "wings" present in the comparison simulation are not 
seen. Instead, the outer envelope material is stripped 
away, allowing the remaining compressed core to ex- 
pand adiabatically into the surrounding medium. Also, 
since we run these simulations for more cloud-crushing 
times, we see the continued destruction imparted upon 
the cloud by the constant inflow of material as it gets 
progressively more shredded. After approximately five 
cloud-crushing times, the relative velocities between the 
cloud material and shock material drop to ~0.1wshock as 
the cloud becomes entrained in the flow. By t ^ St^c, the 
cloud has become well mixed with the inflowing material, 
with density variations on the few percent level. 



8 



Silvia, Smith, & Shull 






0.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 



1 T" 



~1 T" 




"^0 

Fig. 5. — Density (top) and temperature (bottom) profiles for ID (solid line), 2D (dashed line), and 3D (dotted line) comparison 
simulations. The x-axis is presented in code units. From left to right, the panels correspond to times t = 0.23tcc, t = 0.72tcc, t = 1.02tcc, 
and t = 2 Altec- 



One of the most striking differences between the sim- 
ulations, with and without cooling, is the formation 
of cold, dense nodules in the simulations with cooling 
turned on. The nodules appear to be amplifications 
of over-densities created in the cloud. The regions of 
high density have sufficiently shortened cooling times and 
therefore cool and collapse approximately isothermally. 
As a result, the dense knots are able to persist through 
their bombardment with inflowing material for a longer 
fraction of the simulation duration. 

Figure |6] shows the evolution of the ejecta clump to 
roughly one-third of the full simulation time. Beyond 
that point, the cloud continues to expand and mix with 
the shocked medium. The formation of high density nod- 



ules is apparent in the r ight column of Figure [6j We note 
that Raga et al. (20071 found similar dense protrusions 
in their iJU simulations of radiative, interstellar bullet 
flow. The number of nodules that form tends to result 
from the balance between shock-heating and cooling. For 
higher shock velocities, and hence higher temperatures, 
the number of persisting over-densities decreases. 

When we run the cooling simulations with an addi- 
tional level of refinement, we generally find these dense 
knots to be smaller and higher in number. This result is 
expected, as the turbulence created by the shock pass- 
ing through the cloud will naturally cascade down to the 
smallest physical scale, which is determined by the reso- 
lution of the simulation. Additionally, these nodules are 
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TABLE 1 

Simulations Parameters 



Simulation 


Cloud Envelope ^ 


X 


^shock 


M 


^CC 


Cooling 








(km s^-'-) 




(yrs) 




1 


(jaussian 


100 


1000 


27.83 


31.69 


None 


2 


Gaussian 


100 


1000 


27.83 


31.69 


SW87 


3 


Gaussian 


100 


3000 


83.49 


10.56 


None 


4 


Gaussian 


100 


3000 


83.49 


10.56 


SW87 


5 


Gaussian 


100 


5000 


139.15 


6.34 


None 


6 


Gaussian 


100 


5000 


139.15 


6.34 


SW87 


7 


Gaussian 


1000 


1000 


8.80 


100.21 


None 


8 


Gaussian 


1000 


1000 


8.80 


100.21 


SW87 


9 


Gaussian 


1000 


3000 


24.60 


33.40 


None 


10 


Gaussian 


1000 


3000 


24.60 


33.40 


SW87 


11 


Gaussian 


1000 


5000 


44.00 


20.04 


None 


12 


Gaussian 


1000 


5000 


44.00 


20.04 


SW87 


13 


Power-law 


1000 


1000 


8.80 


100.21 


None 



^ For the Gaussian and Power-law envelopes, Tcorc — O.Tr ^i^..^ 

^ For the coohng, "SWS7" represents analytic formula by|Sarazin 8z White| ||l987| 
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Fig. 6. — Density projections in the x-z plane weighted by cloud material, ^, for the case of ^ — 1000 and fshock — 1000 km s-l. The left 
column shows the evolution for a non-cooling simulation, while the right column shows that for a cooling simulation. From top to bottom, 
the time progresses in the following order: t = 1.8tcc, t = 2.4icc, t = S.Otcc- 
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constantly being eroded by inflowing material, and at 
late times their size will be comparable to the grid-cell 
size. However, despite the difference in size and count, 
these small knots evolve similarly in density and temper- 
ature, which results in nearly identical dust destruction. 

4. DATA POST-PROCESSING 

4.1. Modeling Grain Distributions Using Tracer 
Particles 

As mentioned previously, the clouds in all of our sim- 
ulations are populated with a uniform distribution of 
tracer particles. We use 4096 particles in order to ad- 
equately sample the cloud while simultaneously keeping 
the computational requirements of our post-processing 
analysis manageable. Each individual particle is assigned 
a full distribution of grain radii and therefore represents 
a distinct population of dust. The radius distributions 
used in all of our analy sis come from the values calcu- 
lated by Nozawa et al. (2003) for a core-collapse super- 



nova (COyiNj with a progenitor mass of 20 Mq, which 
have been reproduced in Figure [7j Since we assume that 
the number density of dust grains in the cloud should 
scale with the surrounding gas density, we scale the rel- 
ative abundance of grains in the initial distributions by 
the density of the surrounding medium that the tracer 
particles are initially embedded in. Given this assump- 
tion, /(a) = /o(a)(ptp//3max), where a is the grain radius, 
/o(a) is the unsealed abundance, ptp is the density of the 
surrounding medium for a given tracer particle, and Pmax 
is the value for the most dense medium in which a tracer 
particle is initialized. In general, Pmax = Pcorc- To de- 
termine the distribution of grain radii for all the dust 
contained in the simulation, we sum the radius distribu- 
tions for all 4096 tracer particles. 




a [cm] 

Fig. 7. — Reproduced grain radius distributions for the species 
expected in the unmixed ejecta model of a core-collapse sup ernova 
with a progenitor mass of 20 Mq as calculated by Nozaw a et al.| 
( |2003[ l. In situations where data from [Nozawa e t al. (2UU3 ) ended 
prior to a grain radius of 5xl0~* cm, we carried out a linear ex- 
trapolation to smaller grain radii to make the data continuous. 

We can also use the tracer particles and their respective 
grain distributions to calculate a proxy for the amount 
of dust mass in our simulation. First, we define the 
number density of grains between a — da and a to be 



n{a) = f{a)5a, where 5a is set by the number of bins 
used to track the distribution. Second, we convert this 
number density to the mass proxy by multiplying by a^ . 
By summing the amount of mass in each radius bin, we 
can determine the amount of mass contained in a given 
tracer particle and then sum tracer particles to get the 
total dust mass in the simulation. Rather than work in 
absolutes and include the densities of particular grain 
species, we leave this mass proxy as a relative value and 
track how it evolves from t = Xo t — ifinai- 

In the same way that it was difficult to keep all of 
the cloud material inside the finite computational do- 
main, some fraction of the tracer particles also flow off 
the grid, particularly ones that start out in the cloud 
envelope. Since we have no means of tracking the envi- 
ronment that those tracer particles would exist in after 
the point where they leave the bounds of the simula- 
tion, we exclude them from our sputtering calculations 
for the full duration of the run. In some simulations 
the combination of input parameters leads to vigorous 
cloud expansion that cannot be confined to our compu- 
tational domain for the desired duration. This results 
in situations in which large fractions of the tracer parti- 
cles (>50%) flow beyond the box boundaries. To address 
this fact, we halt the post-processing analysis at the point 
when the dust mass contained in the lost tracer particles 
exceeds 5% of the initial dust mass. While we would like 
to carry out our analysis for all simulations on the same 
timescales, creating simulation domains that are large 
enough to contain the cloud's expansion for all of pa- 
rameter space was too computationally expensive to be 
feasible. Despite this limitation, we still find our results 
up to the 5% cut-off to be robust. 

4.2. Dust Sputtering 

Having established our initial grain distributions, we 
post-process the grains using generalized fi fth-order poly- 
nomia l fits to thermal sputtering from [Tielens et al.| 
( 1994 1 and tabulated thermal sputtering rates. These 
sputtering rates have been reproduced in Figure |8] We 
pass the density and temperature of the environment of 
a given tracer particle to one of the sputtering methods 
and calculate the erosion rate, da/dt, for each time step. 
We multiply the erosion rate by the length of the time 
step and subtract that value from the initial grain radius 
such that flnow — loid ~ (da/dt) At. This method allows 
each population of dust (each tracer particle) to evolve in 
direct response to its environment for the duration of the 
simulation. Using the method mentioned above, we then 
compute the pe rcentage of dust mas s lost as a function 



of time. While Tielens et al. (19941 provide sputtering 
rates for SiC an d H2U grains, we do not track them our 
analysis because Nozawa et al. ( 2003 ) do not incl ude size 
distributions for these grains. In addition, since Tielens 



et al. (19941 used Si02 as the basis for the sputtering 
rates 01 silicates, we use the grain distribution of Si02 as 
our starting point for silicate sputtering. 

5. RESULTS 

We present the sputtering results for all 13 simulations 
listed in Table [T] Specifically, we study the changes in 
grain radius distributions from their starting to ending 
values and the evolution of the dust mass contained in 
the tracer particles as a function of time. For a condensed 
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Fig. 8. — Our adopted thermal spu ttering rates based on the 
polynomial fits of |Tielens et al.] | |1994[ | (top ) and for an elemental 
composition of Z = Zq as calculated by [Nozawa et al.| | |2006[ l 
(bottom). 

view of the key focus of this paper, dust survival rate, for 
all combinations of simulations and sputtering methods, 
we refer the reader to Table [2] One caveat is that Simu- 
lations 8 (x = 1000, V shock — 1000 km s^^, cooling on), 

9 (x = 1000, Wshock = 3000 km s^^, cooling off), and 10 
(x = 1000, V shock = 3000 km s~^, cooling on) do not 
span the full range of cloud-crushing times because they 
ended up losing more than 5% of their dust mass due 
to tracer particles that left the computational domain. 
However, since the dust destruction in Simulations 9 and 

10 hits a plateau before the simulation analysis ends, the 
results of these simulations can still be fairly accurately 
compared. Unfortunately, Simulation 8 terminates well 
before any such plateau, leading to misleading final dust 
mass values. In addition, we do not devote any significant 
discussion to Simulation 13 (power-law density envelope, 
X = 1000, -Wshock = 1000 km s~^, cooling off) because 
the dust mass evolution proved to be nearly identical to 
Simulation 7 (gaussian density envelope with the same 
initial conditions), as can be seen in Tabled 

5.1. Grain Distributions 

In order to get a handle on the nature of the dust 
grain sputtering, we first inspect the changes in the dis- 



tributions of grain radius. We find that, regardless of 
species, grains with radii less than a few times 10~® cm 
are sputtered to much lower radii and are often com- 
pletely destroyed. However, given the time scale of our 
simulations, larger grains are capable of surviving the full 
duration of the sputtering, but are shifted to lower radii. 
For simplicity, we defined destroyed grains as those that 
are sputtered to less than 5 x 10~® cm. 

In the figures that follow we look at the evolution of 
grain radius distribution for various locations in param- 
eter space and for two of the three sputtering rate op- 
tions. In each panel, the solid black line illustrates the 
initial grain distribut ions as set by the distributions from 
Nozawa et al. (2003). For the colored lines, solid lines 
correspond to simulations with the slowest shock veloci- 
ties, Wshock — 1000 km s~^, dashed lines are the simula- 
tions with shock velocities of Wshock — 3000 km s^^, and 
dotted lines are the simulations with the highest velocity 
shocks, Ushock = 5000 km s^^. The colors represent our 
different options for cooling, where red highlights runs 
with no cooling and blue highlights runs with th e an- 
alytic cooling function of iSarazin & White (1987). All 



figures adhere to the same scheme. 

In Figures [9| and [T0| we see the evolutio n in three grain 
specie s using the sputtering rates from |Tielens et al.| 
( 1994 1 for density contrasts x = 100 and x = 1000 



respectively. In both the low-density and high-density 
clouds, the higher shock velocities are more effective 
at reducing the size of large grains. In addition, for 
X — 1000, the increased density leads to even more dust 
destruction, with a severe drop off in the small grains as 
a function of shock velocity. When we inspect the grain 
distributions as a function of cooling (red vs. blue), we 
see little difference between the simulations. Since the 
differences are more apparent when looking at dust mass, 
we discuss these differences in the following section. 

We show the resul ts of sputtering using the rates from 
Nozawa et al. (2006) for over-densities of 100 and 100 in 
Figures' TT anciT^ respectively. In the case of x = 100, 
the number of small grains sputtered out of existence 
does not show a clear correlation with shock velocity. In 
some cases, a slower shock destroys more grains than a 
faster shock. This is likely a result of the cloud being 
shredded more rapidly in the high-velocity cases, which 
reduces the density of the cloud material and shuts off 
sputtering. There appears to be no significant depen- 
dence on whether or not cooling was employed in the 
simulation. However, we predict that, with higher metal 
abundances, the differences between cooling and non- 
cooling simulations should become more pronounced. 
For the simulations with x = 1000, grain sputtering de- 
pends strongly on shock velocity, as all grains sizes are 
sputtered more heavily by higher velocity shocks. Again, 
we do not see a clear distinction between cooling and 
non-cooling simulations, even in higher density regimes. 
There does seem to be a large discrepancy for the slow- 
est shock velocity, but this stems from the fact that the 
cooling simulation for that shock velocity did not make 
it to ten cloud-crushing times. The differences between 
simulations with and without cooling can better be seen 
in Section [5.2[ when we look at evolution in dust mass. 

5.2. Dust Mass Evolution 
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Fig. 9. — Grain radius distributions for three of the grain species in |Tielens et al.| ( |1994| l for simulations with over-density x = 100. The 
black line represents the initial distri bution, while the colored lines indicate ftnal distributions. Red signifies simulations with no cooling 
and blue shows the analytic cooling of parazin &; White| |1987jl . Solid colored lines indicate simulations with fshoclc — 

1000 km s-\ dashed 

are v^hnri, = 3000 km s~^, and dotted are Vghodi = 5000 km s~^. The grains in this figure were sputtered using the rates from |Tielens| 
|et al.||l994l l. 




Using the sputtered grain distributions, we convert 
them into relativ e du st masses using the method de- 
scribed in Section |4.1[ In the following figures, we con- 
tinue with the same color and line style convention, but 
instead of radius distributions we consider the fractional 
dust mass of each grain species as a function of cloud- 
crushing times. In these plots, any lines that end before 
10 cloud-crushing times are terminated because more 
than 5% of the original dust mass contained in the parti- 
cles was lost due to particles that crossed the simulation 
boundaries. In addition, for the lower density cloud, ef- 
fectively zero sputtering of iron grains occurred, and we 
have omitted those plots. 

We again look first at the results w hen we use the 
sputtering rates of 'Tielens et al. (1994). Following the 
previous order. Figur e |13| shows tne mass evolution for 
X = 100, and Figure |14| shows the same for x = 1000. 
For the low over-density cloud, with the exception of Fe 
grains, which were negligibly sputtered, we find the scat- 
ter in percentage of dust remaining at the end of the 
simulation is rather small. However, the path that each 
simulation takes to get to its final dust mass varies based 
on the shock velocity. The high- velocity shock runs show 
a rapid drop in dust mass and then dust destruction 
slows considerably, while the dust mass in the low ve- 
locity shock runs gradually falls to its final level. In the 
higher over-density cloud, the dependence on shock ve- 
locity is much more pronounced. Higher shock velocities 



lead to greater dust destruction - as much as 50% addi- 
tional destruction for certain grain species. Comparing 
the red and blue lines in the cooling versus non-cooling 
simulations, we only see differences in the higher den- 
sity cloud, after the initial drop in dust mass. At later 
times, we detect a slight splitting in the dust mass evolu- 
tion tracks in which the simulations with cooling turned 
on appear to destroy a small fraction of additional dust. 
Our explanation for the diverging behavior is that if the 
medium is capable of cooling, it will slightly drop the 
pressure in the expanding cloud and keep the medium 
in which dust is embedded at a higher density. Since 
the cloud has been raised to high shock temperatures at 
late times, where the sputtering rates are relatively con- 
stant with temperature, differences in density will play 

the dominant role in controlling sputtering^ 

We carry out the same inspection for Noz awa et al.| 
p006| in Figure [15] (x = 100) and Figure [l6](x = 1000). 
We irnmediately see that, similar to the results from the 
polynomial sputtering rates, the variation in final dust 
mass is fairly small for the low-density cloud. We also 
see that the trend of faster sputtering at early times 
for higher velocities shocks holds true in these cases as 
well. One particular difference bet ween the results from 
these sputtering rates and those of Tielens et al. ( |1994 1 
is that, for the high- velocity shocks, the fractional dust 
mass seems to plateau at late times, as where in the ear- 
lier figures the dust mass evolution still has a negative 
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Fig. 11. — Grain radius distributions for the six of the grain species in jNozawa et al.] ( |2003| l for simulations with x = 100. The line style s 
and colors have the same meaning as Figure [o] The grains in this figure were sputtered usmg the Z = Zq rates from|Nozawa et al.| ||2006ll. 
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Fig. 12. — Same as Figure [TT] but for over-density x — 1000. 
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Fig. 13. — Dust mass evolution for two of the grain species in |Tielens et al.| ( |1994^ fo r simulations with v = 100. Colors and lines styles 
are the same as Figure [g] The grains in this figure were sputtered using the rates lrom |Tielens et al.| ( |T994| I. 
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Fig. 14. — Same as Figure [13] but for over-density x = 1000. In addition, we include the evolution of Fe grains. 
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Fig. 15. — Dust mas s ev olution for five of the grain species in |Nozawa et al7| ||2003| for simulat ions with y = 100. C olors and line styles 
are the same as Figure [13] The grains in this figure were sputtered using the Z = Zq rates from [Nozawa et al.| | |2006l l. 
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Fig. 16. — Same as Figure [15] but for over-density x = 1000. In addition, we include the evolution of Fe grains. 



slope at the end of the simulation. This may indicate 
that longer simulation run times are required for the 
simulations that did not stop evolving. For the higher 
density density cloud, we again see that the final dust 
mass fractions are much more w i dely v aried. For the 



sputtering rates of Nozawa et al. ( 2006 1 , we see a non- 



negligible amount ot iron sputtering tor the medium ve 
locity shock simulations of the high-density cloud which 
is not present in the polynomial sputtering rate calcula- 
tions. We also see the same effects of cooling in the higher 
density cloud mentioned previously. In addition, one in- 
teresting difference between the non-cooling and cooling 
simulations shows up most prominently in the evolution 
of AI2O3 (see Figure 16 1. The non-cooling simulation 
appears to have a more rapid fall off in dust mass for 
the case of Wshock = 1000 km s^^ than the cooling sim- 
ulation. However, even though the cooling simulation 
did not reach ten cloud-crushing times, it seems unlikely 
that it would not have also arrived at total AI2O3 grain 
destruction. Since this particular dust species consists 
of only small grains, it is possible that the slight drop 
in plasma temperature caused by the cooling is able to 
slow the grain destruction, despite the increased density 
effect mentioned above. 

Finally, in an effort to summarize the results of all 13 
simulations. Tabic [2] lists the final dust mass fraction for 
all species for all simulations. In reviewing these results, 
wc find that grain species that have initial distributions 
weighted toward large grain radii have higher survival 
rates over 10 cloud-crushing times. Furthermore, since 
the sputtering rates scale linearly with density, simu- 
lations with higher density clouds all show substantial 
sputtering. The higher density clouds also tend to show 
a much stronger scaling with shock velocity than do their 



low-density counterparts. There appears to be no strong 
dependence on the two sputtering rate methods, most 
differences are on the tenths of a percent to few percent 
level. While the simulations including cooling did show 
small deviations from those simulations without cooling 
in the dust mass evolution plots, in most cases the end 
results for surviving dust mass fraction are indistinguish- 
able. 

6. SUMMARY AND FUTURE WORK 

We summarize the above the results into a few key 
points: 

• Owing to linear density dependence in the sputter- 
ing rates, the high density clouds destroy signifi- 
cantly higher percentages of dust. In some cases, 
as much as 70% more dust is destroyed, with total 
destruction fractions as high as 80-100%. 

• The relative velocity between the reverse shock 
and the ejecta clump strongly infiuences the total 
dust destruction. For high-density ejecta clumps, 
the highest velocity shock results in an additional 
~50% mass loss compared to the lowest velocity 
shock for some grain species. 

• The initial size distribution of dust grains has a 
considerable impact on the survival rate. Species 
with appreciable mass in grains above 0.1 /xm have 
the highest survival rate, losing only 30% or less of 
their mass. Grains that start out at less than 0.1 
/im are often completely destroyed. 

• Total dust destruction varies widely across grain 
species and can be high for some portions of pa- 
rameter space. For high-density ejecta clumps and 
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TABLE 2 
Final Dust Mass Fraction 



Sim. 



AI2O3 



C 



Mg2Si04 



1 

2 

3 

4 

5 

6 

7 

8* 

9 

10 

11 

12 

13 



0.034 
0.036 
0.030 
0.029 
0.054 
0.052 
0.000 
0.107 
0.000 
0.001 
0.001 
0.001 
0.000 



0.963 
0.964 
0.949 
0.949 
0.948 
0.947 
0.880 
0.941 
0.771 
0.762 
0.643 
0.644 
0.882 



/ 0.962 
/ 0.963 
/ 0.950 
/ 0.950 
/ 0.967 
/ 0.967 
/ 0.889 
/ 0.933 
/ 0.738 
/ 0.728 
/ 0.617 
/ 0.620 
/ 0.890 



Sim. MgSiOs 



SiOa 



0.888 
0.891 
0.811 
0.811 
0.846 
0.844 
0.637 
0.861 
0.401 
0.408 
0.177 
0.181 
0.637 

MgO 



1 

2 

3 

4 

5 

6 

7 

8* 

9 

10 

11 

12 

13 



0.423 
0.428 
0.334 
0.333 
0.394 
0.392 
0.139 
0.452 
0.044 
0.054 
0.014 
0.015 
0.140 



0.932 / 
0.934 / 
0.875 / 
0.875 / 
0.799 / 
0.798 / 
0.676 / 
0.930 / 
0.593 / 
0.598 / 
0.240 / 
0.237 / 
0.676 / 



0.916 
0.919 
0.839 
0.840 
0.870 
0.869 
0.674 
0.893 
0.378 
0.453 
0.197 
0.202 
0.675 



0.838 
0.842 
0.774 
0.775 
0.811 
0.809 
0.615 
0.813 
0.445 
0.382 
0.181 
0.185 
0.615 



Sim. 


Si 


FeS 


Fe 


1 


0.998 


0.626 


1.000 / 1.000 


2 


0.998 


0.631 


1.000 / 1.000 


3 


0.998 


0.520 


1.000 / 1.000 


4 


0.998 


0.519 


1.000 / 1.000 


5 


0.998 


0.540 


1.000 / 1.000 


6 


0.998 


0.537 


1.000 / 1.000 


7 


0.997 


0.314 


1.000 / 0.999 


8* 


0.996 


0.674 


1.000 / 0.998 


9 


0.971 


0.171 


0.999 / 0.962 


10 


0.965 


0.198 


0.997 / 0.961 


11 


0.853 


0.046 


0.932 / 0.789 


12 


0.854 


0.047 


0.928 / 0.790 


13 


0.998 


0.315 


1.000 / 0.999 



Note. — For species with t wo values, they show 
resu lts from sputtering r ates of |Tielens et aT!] ( |1994j | 
and |NQzawa et al.| | [2006[ l, res pectively. Specie s with 
only one value show just the [Nozawa et al. ] l|2006|l 
sputte ring rates. In the context ot the |TTeTenset''aL] 
^1994[ sputtering, Si02 represents "silicates" . 

The high final dust mass values in Simulation 8 are 
a result of the analysis being truncated prior to ten 
full cloud-crushing times, as mentioned in the text. 

high-velocity shocks, the often studied species of 
C, SiOa, and Fe show 62%, 20%, and 80% survival, 
respectively. 

From an observational perspective, these results may 
have interesting implications for the IR emission seen in 
supernova remnants, which is usually attributed to dust 
that has been form e d and processed by the supernova ex- 
plosion. Rho et al. ( 2008 ) modeled detailed mid-IR spec- 
tra of Cas A to determme the populations of dust that 
are likely present in the remnant. They find a reasonably 
high abundance of Si dust, which agrees with the high 
survival rate of those grains. Additionally, 
et al.| ( 2009 ) analyzed mid-IR spectra for S^ 



Sandstrom| 

R 1E0102- 

7219 and found both a lack of AI2O3 grains as well as a 



large population of amorphous carbon grains with radii 
of 0.1 /xm, in good agreement with our results. Such ob- 
servations are excellent indicators that the results from 
our simulations could be used to predict emission fea- 
tures in future remnant observations. 

Since the plasma in these shock-cloud interactions will 
become enriched with high-mass ions as a result of dust- 
grain sputtering, we will investigate, in future work, the 
additional effects that this might have. We have prelim- 
inary findings that indicate that for cooling simulations 
with enhanced metallicity {Z > 10 Zq), the results di- 
verge further from the non-cooling simulations. Specif- 
ically, the simulations with increased cooling owing to 
higher metal abundances tend to show a greater amount 
of dust destruction. This likely stems from the previ- 
ously mentioned effect in which the shorter cooling time 
leads to the formation and collapse of localized cloud 
over-densities that persist for a great fraction of the sim- 
ulation time. We also predict that, with accurate sput- 
tering rates for enhanced metal abundances, we should 
see greater dust destruction by high mass ions, since the 
sputtering yields for these ions should will be factors of 
10-100 higher. 

While our current method of post-processing the tracer 
particles is sufficient as a first step in tracking dust popu- 
lations in an ejecta clump, we would like to employ more 
active means of following dust erosion. Therefore, to im- 
prove upon current simulations, we plan to sputter the 
dust grains "on-the-fly" so that dust destruction occurs 
during simulation run time. This should provide some- 
thing closer to an instantaneous sputtering rate rather 
than the current low time resolution da/dt, which is ap- 
plied only at each discrete data dump. In addition, it 
will allow for a means of accounting for feedback effects 
that dust sputtering may have on the system. As the 
dust grains are eroded, they reintroduce metals into the 
gas phase which would allow them to act as additional 
coolants and ions for grain sputtering. 

Computing the feedback of metals by grain sputter- 
ing also motivates us to improve our cooling method by 
tracking the exact abundances of a variety of metals and 
then using Cloudy to get cooling rates that directly scale 
with the metal abundance pattern. If we accurately track 
the abundances of the metals, we can increase the metal 
abundances in the hot plasma as the grains are sput- 
tered. Providing the link between grain sputtering and 
gas cooling could lead to a possible equilibrium between 
grain sputtering and plasma cooling. At this point, col- 
lapsing gas drives up the sputtering rate, which leads to 
increased metal abundances that cool the gas even fur- 
ther, perhaps to low enough temperatures to shut off the 
thermal sputtering. 

Finally, given the nature of the tracer particles used 
to track dust populations, we have only been able to 
model the dust as being directly coupled to the cloud 
material. Therefore, we have only applied thermal sput- 
tering rates and neglect non-thermal sputtering by grain 
motion through the plasma. Dust coupling to the cloud 
may not be an unreasonable assumption if the grains 
are charged and tied to magnetic fields lines that might 
be present in the ejecta material. It would be ideal to 
decouple the grains from the fluid and account for the ef- 
fects of non-thermal sputtering as well as possible grain- 
grain collisions. Additionally, if magnetic field lines are 
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present, it would be im portant to investi gate the effects 
of betatron acceleration ( ShuU 1977 1978 ) as it might en- 
hance sputtering of large grains. Methods for including 
these factors, as well as those mentioned in the preced- 
ing paragraphs, are being considered for inclusion in a 
second paper on this topic. 
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